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Any model of tides is based on a specific hypothesis of how lagging depends on the tidal- 
flexure frequency x- ^or example, Gerstenkorn (1955), MacDonald (1964), and Kaula 
(1964) assumed constancy of the geometric lag angle 5 , while Singer (1968) and Mignard 
(1979, 1980) asserted constancy of the time lag At . Thus, each of these two models was 
based on a certain law of scaling of the geometric lag: the Gerstenkorn-MacDonald-Kaula 
theory implied that 6 ^ i while the Singer- Mignard theory postulated 5 ~ x^- 

The actual dependence of the geometric lag on the frequency is more complicated 
and is determined by the rheology of the planet. Besides, each particular functional 
form of this dependence will unambiguously fix the appropriate form of the frequency 
dependence of the tidal quality factor, Q{x)- Since at present we know the shape of the 
function Q{x) i can reverse our line of reasoning and single out the appropriate actual 
frequency-dependence of the lag, S{x) ■ as within the frequency range of our concern, 
(5 ~ x° ) = 0.2 — 0.4, then 5 ~ x~"- This dependence turns out to be different 
from those employed hitherto, and it entails considerable alterations in the time scales of 
the tide-generated dynamical evolution. Phobos' fall on Mars is an example we consider. 

1 Introduction 

If a satellite is located at a planetocentric position r, it generates a tidal bulge that either 
advances or retards the satellite motion, depending on the interrelation between the planetary 
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Figure 1: A planet and a tide-raising moon. This picture illustrates the case of a satellite located 
below the synchronous orbit, so that its mean motion n exceeds the planet's spin rate ujp , and the 
tidal bulge is lagging. The angular lag defined as 5 = \f\/r = \up x r — v\ will, generally, 
differ from the absolute value of the angle 6i subtended at the planet's centre between the directions 
to the satellite and the bulge. Since in our study we consider an example with a small eccentricity 
and inclination, we make no distinction between 5 and \ . 

spin rate Up and the tangential part of satellite's velocity v divided by r = |r|. It is convenient 
to imagine (as on Fig. [1]) that the bulge emerges beneath a fictitious satellite located at 

rf = r+f , (1) 



where the position lag / is given by 

f = At {L3pxr - v) . (2) 

At is the time lag between the real and fictitious tide-generating satellites, and the inclination 
and eccentricity of the satellite are assumed sufficiently small. 

The fictitious satellite is merely a way of illustrating the time lag between the tide-raising 
potential and the distortion of the body. This concept implies no new physics, and is but a 
convenient metaphor employed to convey that at each instance of time the dynamical tide is 
modeled with a static tide where all the time-dependent variables are shifted back by At , 
i.e., (a) the moon is rotated back by vAt, and (b) the attitude of the planet is rotated 
back by C3p At . From the viewpoint of a planet-based observer, this means that a dynamical 
response to a satellite located at r is modeled with a static response to a satellite located at 
Vf = r — At {v ~ ujp X r) . 
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In this paper, we intend to dwell on geophysical issues - the frequency-dependence of the 
attenuation rate and its consequences. Hence, to avoid unnecessary mathematical complica- 
tions, in the subsequent illustrative example we shall restrict ourselves to the simple case of 
a tide-raising satellite on a near-equatorial near-circular orbit. In this approximation , the 
velocity of the satellite relative to the surface is 

\ijjp X r — v\ = r \ ujp — n\ , (3) 

the principal tidal frequency is 

X = 2\ujp - n\ , (4) 



and the angular lag is 



6 = — \u;pX r - v\ = — X , (5) 
r 2 



n being the satellite's mean motion, and ujp being the planet's spin rate. The factor of two 
emerges in (jlj) since the moon causes two elevations on the opposite sides of the planet. It will 
also be assumed that x^t ^ 1 ? fo^' which reason we shall neglect the second-order difference 
between the expression ([5]) and the angle subtended at the planet's centre between the moon 
and the tidal bulge (rigorously speaking, the sine of the subtended angle is equal to / x r/r^ ). 

The starting point of all tidal models is that each elementary volume of the planet is subject 
to a tide-raising potential, which in general is not periodic but can be expanded into a sum of 
periodic terms. Within the linear approximation introduced by Love, the tidal perturbations of 
the potential yield linear response of the shape and linear variations of the stress. In extension 
of the linearity approximation, it is always implied that the overall dissipation inside the planet 
may be represented as a sum of attenuation rates corresponding to each periodic disturbance: 



where, at each frequency Xi 



( E(\ ) ) E (y.) 

(...) standing for averaging over flexure cycle, E{xJ denoting the energy of deformation 
at the frequency y- , and Q{Xi) being the quahty factor of the material at this frequency. 
Introduced empirically as a means to figleaf our lack of knowledge of the attenuation process 
in its full complexity, the notion of Q has proven to be practical due to its smooth and 
universal dependence upon the frequency and temperature. An alternative to employment of 
the empirical Q factors would be comprehensive modeling of dissipation using a solution of the 
equations of motion, given a rheological description of the mantle (Mitrovica & Peltier 1992; 
Hanyk, Matyska & Yuen 1998, 2000; Moore & Schubert 2000). Though providing a valuable 
yield for a geophysicist, this comprehensive approach may be avoided in astronomy, where only 
the final outcome, the frequency dependence of Q , is important. Fortunately, this dependence 
is already available from observations. 

In this paper we shall restrict ourselves to the simple case of an equatorial or near-equatorial 
satellite describing a circular or near-circular orbit. Under these circumstances only the prin- 
cipal tidal frequency (jlj) will matter. 
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2 The quality factor Q and the geometric lag angle 5. 



During tidal flexure, the energy attenuation through friction is, as ever, accompanied by a phase 
shift between the action and the response. The tidal quality factor is interconnected with the 
phase lag e and the angular lag 5 via 

= tane = tan 25 (8) 

or, for small lag angles, 

^ e = 26 . (9) 

The doubling of the lag is a nontrivial issue. Many authors erroneously state that is 
equal simply to the tangent of the lag, with the factor of two omitted. For example, Rainey & 
Aharonson (2006) assume that Q^^ is equal to the tangent of the geometric lag. As a result, 
they arrive at a value of Q that is about twice larger than those obtained by the other teams. 
In Bills et al. (2005), one letter, 7 , is used to denote two different angles. Prior to equation (24) 
in that paper, 7 signifies the geometric lag (in our notations, 61 ) . Further, in their equations 
(24) and (25), Bills et al. employ the notation 7 to denote the phase lag (in our notations, e , 
which happens to be equal to 26i ). With this crucial caveat. Bills' equation Q = 1/ tanj is 
correct. This inaccuracy in notations has not prevented Bills et al. (2005) from arriving to a 
reasonable value of the Martian quality factor, 85.58 ± 0.37. (A more recent study by Lainey 
et al. (2007) has given a comparable value of 79.91 ± 0.69 .) 

In the Appendix, we offer a simple illustrative calculation, which explains whence this factor 
of two stems. 

Formulae ([Hl-E]) look reasonable: the higher the quality factor, the lower the damping rate 
and, accordingly, the smaller the lag. What looks very far from being OK are the frequency 
dependencies ensuing from the assertions of S being either constant or linear in frequency: 
the approach taken by Gerstenkorn (1955), MacDonald (1964), and Kaula (1964) implies that 
Q ~ X° , while the theory of Singer (1968) and Mignard (1979, 1980) yields Q ~ > neither 
option being in agreement with the geophysical data. 

3 Dissipation in the mantle. 

3.1 Generalities 

Back in the 60s and 70s of the past century, when the science of low-frequency seismological 
measurements was yet under development, it was widely thought that at long time scales 
the quality factor of the mantle is proportional to the inverse of the frequency. This fallacy 
proliferated into planetary astronomy where it was received most warmly, because the law 
Q ~ 1/x turned out to be the only model for which the linear decomposition of the tide 
gives a set of bulges displaced from the direction to the satellite by the same angle. Any 
other frequency dependence Q{x) entails superposition of bulges corresponding to the separate 
frequencies, each bulge being displaced by its own angle. This is the reason why the scaling 
law Q ~ 1/x, long disproved and abandoned in geophysics (at least, for the frequency band 
of our concern), still remains a pet model in celestial mechanics of the Solar system. 

Over the past twenty years, considerable progress has been achieved in the low-frequency 
seismological measurements, both in the lab and in the field. Due to an impressive collective 
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effort undertaken by several teams, it is now a firmly established fact that for frequencies down 
to about ~ 1 yr~^ the quality factor of the mantle is proportional to the frequency to the power 
of a positive fraction a . This dependence holds for all rocks within a remarkably broad band 
of frequencies: from several MHz down to about 1 yr^^ . 

At timescales longer than 1 yr, all the way to the Maxwell time (about 100 yr), attenuation 
in the mantle is defined by viscosity, so that the quality factor is, for all minerals, well approx- 
imated with r]x/M , where f] and M are the shear viscosity and the shear elastic modulus of 
the mineral. Although the values of both the viscosity coefficients and elastic moduli greatly 
vary for different minerals and are sensitive to the temperature, the overall quality factor of 
the mantle at such long timescales still is linear in frequency. 

At present there is no consensus in the seismological community in regard to the time scales 
exceeding the Maxwell time. One viewpoint (incompatible with the Maxwell model) is that 
the linear law Q ~ x extends all the way down to the zero- frequency limit (Karato 2007). An 
alternative point of view (prompted by the Maxwell model) is that at scales longer than the 
Maxwell time we return to the inverse-frequency law Q ~ 1/x . 

All in all, we have: 

For 10^ Hz > X > ly^"^ : Q ~ , with a = 0.2 - 0.4 ( 0.2 for partial melts) . (10) 
For 1 yr'^ > X > 10"^ yr'^ : Q ~ X • (H) 

For 10^^ yr^^ > x '■ arguably, it is still Q ~ x • ( Or maybe Q ~ 1/x ?) (12) 

Fortunately, in practical calculations of tides in planets one never has to transcend the Maxwell 
time scales, so the controversy remaining in (fT2|) bears no relevance to our subject. We leave for 
a future study the case of synchronous satellites, the unique case of the Pluto-Charon resonance, 
or the binary asteroids locked in the same resonance. Thus we shall avoid also the frequency 
band addressed in flTTl) . but shall be interested solely in the frequency range described in (ITU]) . 
It is important to emphasise that the positive-power scaling law fllOj) is well proven not only for 
samples in the lab but also for vast seismological basins and, therefore, is universal. Hence, this 
law may be extended to the tidal friction - validity of this extension will be discussed below in 
subsection 3.4.1 

Below we provide an extremely condensed review of the published data whence the scaling 
law flTU]) was derived by the geophysicists. The list of sources will be incomplete, but a full 
picture can be obtained through the further references contained in the works to be quoted 
below. For a detailed treatment, see Chapter 11 of the book by Karato (2007) that contains a 
systematic introduction into the theory of and experiments on attenuation in the mantle. 

3.2 Circumstantial evidence: attenuation in minerals. 
Laboratory measurements and some theory 

Even before the subtleties of solid-state mechanics with or without melt are brought up, the 
positive sign of the power a in the dependence Q ~ x" ^i^-Y be anticipated on qualitative 
physical grounds. For a damped oscillator obeying z + 2Pz + x'^^ = , the quality factor 
is equal to x/(2/5) , i.e., Q ~ X • 
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Solid-state phenomena causing attenuation in the mantle may be divided into three groups: 
the point-defect mechanisms, the dislocation mechanisms, and the grain-boundary ones. 

Among the point-defect mechanisms, most important is the transient diffusional creep, i.e., 
plastic flow of vacancies, and therefore of atoms, from one grain boundary to another. The flow 
is called into being by the fact that vacancies (as well as the other point defects) have different 
energies at grain boundaries of different orientation relative to the applied shear stress. This 
anelasticity mechanism is wont to obey the power law Q ~ with a ~ 0.5 . 

Anelasticity caused by dislocation mechanisms is governed by the viscosity law Q ~ x 
valid for sufficiently low frequencies (or sufficiently high temperatures), i.e., when the viscous 
motion of dislocations is not restrained by the elastic restoring stress. (At higher frequencies 
or/and lower temperatures, the restoring force "pins" the defects. This leads to the law Q ~ 
(1 + T^x^)''"~^X~^ ! parameter r being the relaxation time whose values considerably vary 
among different mechanisms belonging to this group. As the mantle is warm and viscous, we 
may ignore this caveat.) 

The grain-boundary mechanisms, too, are governed by the law Q ~ x"; though with a lower 
exponent: a ~ 0.2 — 0.3. This behaviour gradually changes to the viscous mode [a = 1) at 
higher temperatures and/or at lower frequencies, i.e., when the elastic restoring stress reduces. 

We see that in all cases the quality factor of minerals should grow with frequency. Accord- 
ingly, laboratory measurements confirm that, within the geophysically interesting band of x ? 
the quality factor behaves as Q ~ with a = 0.2 — 0.4 . Such measurements have been 
described in Karato & Spetzler (1990) and Karato (1998). Similar results were reported in the 
works by the team of I. Jackson - see, for example, the paper (Tan et al. 1997) where numerous 
earlier publications by that group are also mentioned. 

In aggregates with partial melt the frequency dependence of Q keeps the same form, with 
a leaning to 0.2 - see, for example, Fontaine et al. (2005) and references therein. 

3.3 Direct evidence: attenuation in the mantle. 
Measurements on seismological basins 

As we are interested in the attenuation of tides, we should be prepared to face the possible 
existence of mechanisms that may show themselves over very large geological structures but 
not in small samples explored in the lab. No matter whether such mechanisms exist or not, we 
would find it safer to require that the positive-power scaling law Q ~ x" 5 even though well 
proven in the lab, must be propped up by direct seismological evidence gathered over vast zones 
of the mantle. Fortunately, such data are available, and for the frequency range of our interest 
these data conform well with the lab results. The low-frequency measurements, performed by 
different teams over various basins of the Earth's upper mantle, agree on the pivotal fact: the 
seismological quality factor scales as the frequency to the power of a positive fraction a - see, 
for example, Mitchell (1995), Stachnik et al. (2004), Shito et al. (2004), and further references 
given in these sources^ 

^ So far, Figure 11 in Flanagan & Wiens (1998) is the only experimental account we know of, which only 
partially complies with the other teams' results. The figure contains two plots depicting the frequency depen- 
dencies of 1/Qshear and 1/Q compress ■ While the behaviour of both parameters remains conventional down to 
10~^Hz, the shear attenuation surprisingly goes down when the frequency decreases to 10"'^ Hz. Later, one of 
the Authors wrote to us that "Both P and S wave attenuation becomes greater at low frequencies. The trend 
towards lower attenuation at the lowest frequencies in Fig. 11 is not well substantiated.''^ (D. Wiens, private 
communication) Hence, the consensus on (jlOp stays. 
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3.4 Consequences for the tides 



3.4.1 Tidal dissipation vs seismic dissipation 

For terrestrial planets, the frequency-dependence of the Q factor of bodily tides is similar to 
the frequency-dependence (TTUl- fTT]) of the seismo logical Q factor. This premise is based on the 
fact that the tidal attenuation in the mantle is taking place, much like the seismic attenuation, 
mainly due to the mantle's rigidity. This is a nontrivial fact because, in distinction from 
earthquakes, the damping of tides is taking place both due to nonrigidity and self-gravity of 
the planet. Modeling the planet with a homogeneous sphere of density p, rigidity /i , surface 
gravity g, and radius R, Goldreich (1963) managed to separate the nonrigidity-caused and 
self-gravity-caused inputs into the overall tidal attenuation. His expression for the tidal quality 
factor has the form 

Qo being the value that the quality factor would assume were self-gravity absent. To get an idea 
of how significant the self-gravity-produced input could be, let us plug there the mass and radius 
of Mars and the rigidity of the Martian mantle. For the Earth's mantle, /x = 65 -i- 80 GPa. 
Judging by the absence of volcanic activity over the past hundred(s) of millions of years of 
Mars' history, the temperature of the Martian upper mantle is (to say the least) not higher 
than that of the terrestrial one. Therefore we may safely approximate the Martian /i with 
the upper limit for the rigidity of the terrestrial mantle: /x = 10^^ Pa. All in all, the relative 
contribution from self-gravity will look as 

— 11 3 —1 —2 23 2 

2 gpR _ 6 7 M2 _ 1 ( 6.7 X 10 m kg s ) ( 6.4 x 10 kg ) _ -2 
19 ~ir ~ n^JW m (10" Pa) (3.4 X 10' m)' ^'^ 

7 denoting the gravity constant. This, very conservative estimate shows that self-gravitation 
contributes, at most, several percent into the overall count of energy losses due to tides. This 
is the reason why we extend to the tidal Q the frequency-dependence law measured for the 
seismic quality factor. 

3.4.2 Dissipation in the planet vs dissipation in the satellite 

A special situation is tidal relaxation toward the state where one body shows the same side to 
another. Numerous satellites show, up to librations, the same face to their primaries. Among 
the planets, Pluto does this to Charon. Such a complete locking is typical also for binary as- 
teroids. A gradual approach toward the synchronous orbit involves ever- decreasing frequencies, 
eventually exceeding the limits of equation ffTTl) and thus the bounds of the present discussion. 
Mathematically, this situation still may be tackled by means of flTUl) until the tidal frequency 
X decreases to lyr~^, and then by means of flTT]) while x remains above the inverse Maxwell 
time of the planet's material. Whether the latter law can be extended to longer time scales 
remains an open issue of a generic nature that is not related to a specific model of tides or to a 
particular frequency dependence of Q . The generic problem is whether we at all may use the 
concept of the quality factor beyond the Maxwell time, or should instead employ, beginning 
from some low x 1 ^ comprehensive hydrodynamical model. In the current work, we address 
solely the satellite-generated tides on the planet. The input from the planet-caused tides on 



the satellite will be considered elsewhere. The case of Pluto will not be studied here either. 
Nor shall we address binary asteroids. (Since at present most asteroids are presumed loosely 
connected, and since we do not expect the dependencies ( ITO] - [11]) to hold for such aggregates, 
our theory should not, without some alterations, be applied to such binaries.) 

Thus, since we are talking only about dissipation inside the planet, and are not addressing 
the exceptional Pluto-Charon case, we may safely assume the tidal frequency to always exceed 
lyr~^. Thence ffTOl) will render, for a typical satellite: 

g ~ , with a = 0.2 - 0.4 . (15) 

Accordingly, ([9]) will entail: 

5 ~ X"" , with a = 0.2 - 0.4 . (16) 

Another special situation is a satellite crossing a synchronous orbit. At the moment of 
crossing, the principal tidal frequency x = 2 jcUp — n\ vanishes. As ([9]) and ( ITOl) yield Q ~ 
(25) ~^ and Q ~ > then we get S ~ x^" with a positive a. Uncritical employment of 
these formulae will then make one think that at this instant the lag 6 grows infinitely, a clearly 
nonsensical result. The quandary is resolved through the observation that the bulge is lagging 
not only in its position but also in its height, for which reason the dissipation rate remains 
finite (Efroimsky 2007). Since in this paper we shall not consider crossing of or approach to 
synchronous orbits, and since the example we aim at is Phobos, we shall not go deeper into 
this matter here. 



3.5 A thermodynamical aside: the frequency and the temperature 

In the beginning of the preceding subsection we already mentioned that though the tidal Q 
differs from the seismic one, both depend upon the frequency in the same way, because this 
dependence is determined by the same physical mechanisms. This pertains also to the tem- 
perature dependence, which for some fundamental reason combines into one function with the 
frequency dependence. 

As explained, from the basic physical principles, by Karato (2007, 1998), the frequency and 
temperature dependencies of Q are inseparably connected. The quahty factor can, despite its 
frequency dependence, be dimensionless only if it is a function not just of the frequency per 
se but of a dimensionless product of the frequency by the typical time of defect displacement. 
This time exponentially depends on the activation energy A*, whence the resulting function is 

g ~ [X exp{A*/RT)r . (17) 

For most minerals of the upper mantle. A* lies within the limits of 360 — 540 kJ mol~^. For 
example, for dry olivine it is about 520 kJ mol~^. 

Thus, through formulae f|T71) and IQ, the cooling rate of the planet plays a role in the orbital 
evolution of satellites: the lower the temperature, the higher the quality factor and, thereby, the 
smaller the lag 6 . For the sake of a crude estimate, assume that most of the tidal attenuation 
is taking place in some layer, for which an average temperature T and an average activation 
energy A* may be introduced. Then from (fT7|l we have: AQ/Q ^ —aA*AT/RT. For a 
reasonable choice of values a = 0.3 and A* = 5.4 x 10^ J /mol, a drop of the temperature 
from To = 2000 K down by AT = 200 K will result in AQ/Q ^ 1 . So a 10 % decrease of 
the temperature can result in an about 100 % growth of the quality factor. 

Below we shall concentrate on the frequency dependence solely. 
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4 Formulae 



The tidal potential perturbation acting on the tide-raising satellite is 



A, 



^5 



(3 (rj-r)' 



-2 --2 

rj r 



) + 



A, 



iy^6 6 



( 5 {vf ■ r) — 3 r J r ) + . 



:i8) 



where r = |r| and Vf = \rf\, while the constants are given by 

k2GmR^ , _ k^GmW 



A, 



A, 



(19) 



kn being the Love numbers. (For derivation of (fT8|) see, for example, MacDonald (1964) and 
literature cited therein.) 

Three caveats will be appropriate at this point. First, to (ITSl) we should add the potential 
due to the tidal distortion of the moon by the planet. That input contributes mainly to the 
radial component of the tidal force exerted on the moon, and entails a decrease in eccentricity 
and semi-major axis (MacDonald 1964). Here we omit this term, since our goal is to clarify the 
frequency dependence of the lag. Second, we acknowledge that in many realistic situations the 
/cs and sometimes even the ^4 term is relevant (Bills et al. 2005). With intention to keep these 
inputs in our subsequent work, here we shall restrict our consideration to the leading term only. 
Hence the ensuing formula for the tidal force will readj^ 



3 ^2 G 



2,k2Gm^ R^ 



r 

r 



-* 2 
r r 



f 

r 



fr^ - 2r (f-f 



(20) 



The third important caveat is that in our further exploitation of this formula we shall take 

— * 

into account the frequency- dependence of the lag / , but not of the parameter k2 . While the 

— * — * 

dependence f{x) '^iH be derived through the interconnection of / with 5{x) ^^id therefore 
with Q{x) , the value of k2 will be asserted constant. That the latter is acceptable can be 
proven through the following formula obtained by Darwin (1908) under the assumption of the 
planet being a Maxwell body (see also Correia & Laskar 2003): 



1 + ^2/^2 



''fluid 



1 + ix'vVf^') (1 + l9fi/i2gpR)y 



Here kji^-^ is the so-called fluid Love number. This is the value that k2 would have assumed 
had the planet consisted of a perfect fluid with the same mass distribution as the actual planet. 
Notations /i , p , g , and g stand for the rigidity, mean density, surface gravity, and the 
radius of the planet. For these parameters, we shall keep using the estimates from subsection 
3.4. The letter rj signifies the viscosity. Up to an order or two of magnitude, its value may be 



approximated, for a terrestrial planet's mantle, with 10^^ kg/ (m- s 
(x • 10"" 

the frequency-dependence in Darwin's formula will cancel out. 



This will yield: x V / p 



^^s)^ , wherefrom we see that in all realistic situations pertaining to terrestrial planets 



Thus we shall neglect the 



Be mindful that 0(/ /r^) ^ 0{S'^) <C 1 
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frequency-dependence of the Love number ^2 (but shall at the same time take into account 
the frequency-dependence of Q , for it will induce frequency-dependence of all three lags). 
The interconnection between the position, time, and angular lags, 

6^ = At i|d;,xr - t;| =^x , (21) 
can be equivalently rewritten as: 

/ = f r5 = r^Xf , (22) 



where 



cjp X r — v 
ijjp X r — V 



(23) 



is a unit vector pointing in the lag direction. 

Be mindful that we assume the inclination and eccentricity to be small, wherefore the ratio 

Ifl 1 

— = At - \ ujr, X r — V \ 



r 



is simply the tangential angular lag, i.e., the geometric angle subtended at the primary's centre 
between the moon and the bulge. In the general case of a finite inclination or/and eccentricity, 
all our formulae will remain in force, but the lag 6 will no longer have the meaning of the 
subtended angle. 

At this point, it would be convenient to introduce a dimensional integral parameter S 
describing the overall tidal attenuation rate in the planet. The power scaling law mentioned in 
section 3 may be expressed as 

g = , (24) 

where is simply the dimensional factor emerging in the relation Q ~ . As mentioned 
in subsection 3.5, cooling of the planet should become a part of long-term orbital calculations. 
It enters these calculations through evolution of this parameter S . Under the assumption that 
most of the tidal dissipation is taking place in some layer, for which an average temperature T 
and an average activation energy A* may be defined, f|T7|) yields: 



'A* 






~R 







£ = £o exp 

To being the temperature of the layer at some fiducial epoch. The physical meaning of the 
integral parameter £ is transparent: if the planet were assembled of a homogeneous medium, 
with a uniform temperature distribution, and if attenuation in this medium were caused by one 
particular physical mechanism, then £ would be a relaxation time scale associated with this 
mechanism (say, the time of defect displacement). For a realistic planet, £ may be interpreted 
as a relaxation time averaged (in the sense of Q = £°'x°') over the planet's layers and over 
the various damping mechanisms acting within these layers. 

As dlO]) entails 6 ^ 1/(2(5) = (1/2) , then (Ell) necessitates for the position lag: 

/ = r 5 f = ^ r f X"" f , (25) 
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and for the time lag: 

At = , (26) 

S being the planet's integral parameter introduced above, and x being a known function (j4]) 
of the orbital variables. Putting everything together, we arrive at 

/ = i (Ex)-' a "S^i-t = i iSxr a ^t' + 0(e) . (27) 

2 1 + e cos u \ ujp X r — V \ 2 \LiJp x r — v \ 

where 

X=2\LUp-n\ . (28) 



The time lag is, according to fl26|) : 

-{a+l) 

At = S {2 £ \ujp - n\) . (29) 

Formulae fl29|) . fl27|l . and fl20|) are sufficient to both compute the orbit evolution and trace the 
variations of the time lag. 



5 The example of Phobos' fall to Mars 

As an illustrative example, let us consider how the realistic dependence Q{x) alters the life 
time left for Phobos. We shall neglect the fact that Phobos is close to its Roche limit, and may 
be destroyed by tides prior to its fall. We also shall restrict the dynamical interactions between 
Phobos and Mars to a two-body problem disturbed solely with the tides raised by Phobos on 
Mars. Thus we shall omit all the other perturbations, like the Martian non-sphericity and 
precession, or the pull exerted upon Phobos by the Sun, the planets, and Deimos. If, along 
with these simplifications, we assume the eccentricity and inclination to be small, then we shall 
be able to describe the evolution of the semi-major axis by means of the following equation 
(Kaula 1964, p. 677, formula 41): 

da 3 ko Gm , ^ 

dt Q y/G {Mo + m) all/2 ' 

with Mo and m denoting the masses of Mars and Phobos. This equation can be solved 
analytically, provided the quality factor Q is set constant (as in Kaula 1964). The solution is: 

/ ^ \ 2/13 

/ 39 k2R^Gm ,0/2 \ 

a{t) = / t + /2 3^ 

^ 2Q^GiM + m) "J 

tto = o.{t)\t=Q being the initial value of Phobos' semi-major axis. 

Unfortunately, neither our model (wherein Q is given by (l24l) ) nor the Singer-Mignard 
model (with Q scaling as 1/x ) admit such an easy analytical solution. This compels us to 
rely on numerics. In a ( quasi )inertial frame centered at Mars, the equation of motion looks: 

(fr G{Mo + m)r 3 G {Mo + m) m R'' \ -.3 ^ ^ / ^ 1 

^ = 7^ y-^' -2r(r./jJ (32) 
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Naturally, its right-hand side consists of the principal, two-body contribution — G{Mo + ■m)r/r^ 
and the disturbing tidal force given by (!20!) . It should be noted, however, that we did not bring 
in here all the terms from (120!) . Following Mignard (1980), we retain in (!32l) only the per- 
turbing terms dependent on /. The other perturbing term (the first term on the right-hand 
side of f l2U]) ) is missing in fl52]) . because it provides no secular input into the semi- major axis' 
evolution. For a proof of this fact see Appendix A. 3 below. 

Phobos' orbital motion obeys the planetary equations in the Euler-Gauss form. In assump- 
tion of i and e being small, the problem conveniently reduces to one equation: 

1 = ^ . (33) 

at n 

where S is the tidal acceleration given by the second term of fl52]) projected onto the tangential 

^ — * — * 

direction of the satellite motion. This direction is defined by the unit vector {H x r)/\H x r\ , 
where H denotes the angular-momentum vector. Thence the said projection reads: 

3 k2G (Mo + m) m r 2. - ^ ^ , (Hxr) 

^- ^ V o ; [- At {ujp X r - v) + 2rAt {r ■ v)] ^ ' 



r 
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Mo ' ' ' ' " " \Hxi^ 



In Mignard (1980) the appropriate expression is given with a wrong sign. This is likely to be 
a misprint, because the subsequent formulae in his paper are correct. 

— * 

In assumption of i and e being negligibly small, (H x r) ■ (^ujp x r) can be approximated 
with na^Up , whereafter (l3l|l gets simplified to 

3k2R'' G {Mo + m)mAt 
S = in - .p) , (35) 

substitution whereof into ( l33l) entails the following equation to integrate: 

da 6 k2 n m At 



dt Mo a^ 



(n - LUp) . (36) 



Our computational scheme was based on the numerical integrator RA15 offered by Everhart 
(1985). The initial value of a , as well as the values of all the other physical parameters entering 
(l36l) . were borrowed from Lainey et al. (2007). These included an estimate of 0.6644 min for 
the present-day time lag At . 

Four numerical simulations were carried out. One of these implemented the Singer-Mignard 
model with a tidal-frequency-independent At . The other three integrations were performed 
for the realistic frequency-dependence fl2^ . with a = 0.2 , 0.3 and 0.4. 

To find the integral parameter S emerging in (12^ . we used the present-day values of At 
and n . The resulting values of S were found to be 1201 x 10^ day rad~^ , 81028 day rad~^ 
and 2104 day rad~^ for a = 0.2 , 0.3 and 0.4 , respectively. We did not take into account the 
strong temperature-dependence of S , leaving this interesting topic for discussion elsewhere. 

Simultaneous numerical integration of equations fl2^ and fl5^ results in plots presented 
on Figures [2] and El The first of these pictures shows the evolution of Phobos' semi-major 
axis from its present value until the satellite crashes on Mars, having descended about 6000 
km. The leftmost curve reproduces the known result that, according to the Singer-Mignard 
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model with a constant At , Phobos should fall on Mars in about 29 MyrEl The next curve 
was obtained not numerically but analytically. It depicts the analytical solution (13T|) available 
for the Gerstenkorn-MacDonald-Kaula model with a constant Q , and demonstrates that this 
model promises to Phobos a longer age, 38 Myr. The three curves on the right were obtained 
by numerical integration of fl29l) and (15^ . They correspond to the realistic rheology with a 
equal to 0.2 , 0.3 and 0.4 . It can be seen that within the realistic model Phobos is expected to 
survive for about 40 — 43 Myr, dependent upon the actual value of a of the Martian mantle. 
This is about 15 Myr longer than within the Singer-Mignard model widely accepted hitherto. 

The difference between the three scenarios shown on Figure [2] stems from the different 
rate of evolution of the lag At in the three theories addressed. Within the Singer-Mignard 
formalism. At stays unchanged through the descent. As can be seen from formula ( l29l) . this 
is equivalent to setting a = — 1 , an assertion not supported by geophysical data. Within the 
Gerstenkorn-MacDonald-Kaula model, the time lag is subject to a gradual decrease described 
by the formula 

At = (37) 

2 \n — ujpl 

under the assumption that Q is constant and is equal to its present-day value Q = 79.91 
determined in (Lainey et al. 2007). Comparison of fl37j) with fj29|) reminds us of the simple fact 
that, in terms of our model, Gerstenkorn-MacDonald-Kaula's theory corresponds to the choice 
of a = , a choice which is closer to the realistic rheology than the Singer-Mignard model. 

In the reahstic model, a is positive and assumes a value of about 0.2 — 0.4 . As a result, the 
time lag is gradually decreasing. However this decrease looks different from that in Gerstenkorn- 
MacDonald-Kaula's model - for their comparison see Figure [31 



6 Conclusions 

As the tidal angular lag 6 is inversely proportional to the tidal Q factor, the actual frequency- 
dependence of both 6 and At is unambiguously defined by the frequency-dependence of 
Q . While in the Gerstenkorn-MacDonald-Kaula theory of tides the geometric lag is assumed 
frequency-independent, in the Singer-Mignard theory it is the time lag that is spared of fre- 
quency dependence. However, neither of these two choices conform to the geophysical data. 

We introduce a realistic tidal model, which permits the quality factor and, therefore, both 
the angular lag 6 and the time lag At to depend on the tidal frequency x ■ The quality 
factor is wont, according to numerous studies, to obey the law Q ~ ? where a lies within 
0.2 — 0.4 . This makes the time lag At not a constant but a function ( l26l) of the principal 
tidal frequency and, through ( l29l) . of the orbital elements of the satellite. The same pertains 
to the angular lag 6. 

Using these tidal-frequency dependencies for the time and angular lags, along with the 
recently updated values of the Martian parameters, we explored the future of Phobos, taking 
into account only the tides raised by Phobos on Mars, but not those caused by Mars on Phobos. 
Our integration shows that Phobos will fall on Mars in 40 — 43 Myr from now. It is up to 50 % 

^ In the paragraph after his formula (18), Mignard (1981) states that "Phobos will end its life in about 
36 million years". Mignard arrived to that number by using an old estimate of 20 deg/cyr^ for the initial 
tidal acceleration. Later studies, like for example Jacobson et al. (1989) and Lainey et al. (2007) have shown 
that this value should be increased to 25.4 deg/cyr^. It is for this reason that our simulation based on the 
Singer-Mignard model gives not 36 but only 29 Myr for Phobos' remaining lifetime. 
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Time (years) 



Figure 2: Evolution of Phobos' semi-major axis, as predicted by different models. The lines (from 
left to right) correspond to the Singer-Mignard model, to the Gerstenkorn-McDonald-Kaula model, 
and to the realistic rheology with a = 0.2 , with a = 0.3 , and with a = 0.4 , correspondingly. 
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Figure 3: Evolution of At computed for three different models. The horizontal green line corresponds 
to the Singer-Mignard model wherein At is set to be constant. The red line corresponds to the 
Gerstenkorn-MacDonald-Kaula model. The dark-blue, violet, and light-blue lines correspond to the 
realistic rheological models with a = 0.2 , with a = 0.3 . and with a = 0.4 , correspondingly. 
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longer than the estimate stemming from the Singer-Mignard model employed in the past. This 
demonstrates that the currently accepted time scales of dynamical evolution, deduced from old 
tidal models, should be reexamined using the actual frequency dependence of the lags. 
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Appendix. 

The goal of this Appendix is threefold. First, we remind the reader why in the first approxima- 
tion the quality factor is inversely proportional to the phase lag. Second, we explain why the 
phase lag is twice the geometric lag angle, as in formulae (I8]-|9]) above. While a comprehen- 
sive mathematical derivation of this fact can be found elsewhere (see the unnumbered formula 
between equations (29) and (30) on p. 673 in Kaula 1964), here we illustrate this counterintu- 
itive result by using the simplest setting. Third, we justify our neglect of the first term in fl2U]) . 



A.l. The case of a near-circular near-equatorial orbit. 



Consider the simple case of an equatorial moon on a circular orbit. At each point of the 
planet, the tidal potential produced by this moon will read 

W = Wo cosxt , (38) 

the tidal frequency being given by 

X = 2 \n — Ljp\ . (39) 

Let g denote the free-fall acceleration. An element of the planet's volume lying beneath the 
satellite's trajectory will then experience a vertical elevation of 

C = — cos(xt - 26) . (40) 

Accordingly, the vertical velocity of this element of the planet's volume will amount to 

u = ( = — X — - sin(X'^ — 2 5) = — x —— (sinx^ cos 25 — cos^t sin 25) . (41) 
g g 

The expression for the velocity has such a simple form because in this case the instantaneous 
frequency x is constant. The satellite generates two bulges - on the facing and opposite sides 
of the planet - so each point of the surface is uplifted twice through a cycle. This entails the 
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factor of two in the expressions fl3^ for the frequency. The phase in fHU]) . too, is doubled, 
though the necessity of this is less evident 

The energy dissipated over a time cycle T = 27r/x, per unit mass, will, in neglect of 
horizontal displacements, be 



AE 



cycle 



^ ^ dW\ f Wo\ dWo 

u I — \ dt = — ( — X — ) — — / cosx^ (sinxt cos 25 — cosx^ dt 



dr 



g / dr 



t=o 



= — X — - ° sin 26 — f cos^ xt d{xt) = " vr sin 26 

g or X Jxt = o g 9r 

while the peak energy stored in the system during the cycle will read: 



(42) 



E 



peak 



T/4 



U 



dW 
dr 



dt 



dr 



cos xt (sin xt cos 26 — cos xt sin 26) dt 



2 a 



g dr 



cos 25 /•^* = -/2 sin 25 /•x* = -/2 

/ cosxt smxt d[xt) / 

X Jxt = o X Jxt = o 



cos^Xt d{xt) 



dW^ 
g dr 



1 vr 

- cos 25 — — sin 25 

2 4 



whence 



Let a; signify a position along the equatorial circumference of the planet. In the absence of lag, the radial 
elevation at a point x would be: 

Wo 



C 



cosfc(a; — V t) , v = Ra 



V being the velocity of the satellite's projection on the ground, R being the planet's radius, and a being simply 
|n — Wpl because we are dealing with a circular equatorial orbit. The value of k must satisfy 

k V — 2 a , i.e., k v = x j 

to make sure that at each x the ground elevates twice per an orbital cycle. The above two formulae yield: 

kR = 2 . 

In the presence of lag, all above stays in force, except that the formula for radial elevation will read: 

C — cosk(x — V t + D) , where D — RS , 
g 

D being the linear lag, and S being the angular one. Since kv = 2, we get: 

cos [k {x — V t + R S-^)] — cos [k x — k v t + k R 5] — cos [k x — {k v t — 2 5)] , 
so that, at some fixed point (say, at a; = 0) the elevation becomes: 

at) = ^ cos(fc vt - 25) . 
g 

We see that, while the geometric lag is 5 , the phase lag is double thereof. 



(43) 
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Q ^ = = ^ tan 25 . (44) 

27r^^,,, 27r i cos25 - | sin25 



The above formulae were written down in neglect of horizontal displacements, approximation 
justified below in the language of continuum mechanics. 

A. 2. On the validity of our neglect of the horizontal displacements 

In our above derivation of the interrelation between Q and 6 , we greatly simplified the 
situation, taking into account only the vertical displacement of the planetary surface, in response 
to the satellite's pull. Here we shall demonstrate that this approximation is legitimate, at least 
in the case when the planet is modeled with an incompressible and homogeneous medium. 

As a starting point, recall that the tidal attenuation rate within a tidally distorted planet 
is well approximated with the work performed on it by the satellit^: 

E = - j pV VW £x (45) 

where p , V and W are the density, velocity, and tidal potential inside the planet. To simplify 
this expression, we shall employ the equality 

pVVW = V ■ (pVW) - WV-Vp - WV-{pV) = V-{pV VW) - W V -Vp + W . (46 

— * 

For a homogeneous and incompressible primary, both the V WVp and dp/dt terms are nil, 
wherefrom 

E=-jpWV-nd^x , (47) 

n being the outward normal to the surface of the planet. We immediately see that, within the 
hydrodynamical model, it is only the radial elevation rate that matters. 

Now write the potential as W = Wo cos(xt) . Since the response is delayed by At, the 

— * 

surface-inequality rate will evolve as V ■ n ~ sin [x{t — ^0 ] • the rest will then be as 
in subsection A.l above. 

A. 3. On the validity of our neglect of the nondissipative tidal potential 

The right-hand side of equation (l32l) consists of the principal part, — G {Mo + m)r/r^ , 
and tidal perturbation terms. These are the second and third terms from the right-hand side 
of fl20p . terms that bear a dependence on / and, therefore, on At. The first term from the 
right-hand side of fl20p lacks such a dependence and, therefore, is omitted in (1321) . The term 
was dropped because it would provide no secular input into the history of the semi-major axis. 
Here we shall provide a proof of this statement. 

The omitted term corresponds to a potential (Mignard 1980): 

k2{Mo + m)GmR^ 

- — 2m:7^ — • ^^^^ 



A small share of this work is being spent for decelerating the planet rotation 
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From the physical standpoint, Uq models the effect of the tidal bulges, assuming their direc- 
tion to coincide with that toward the tide-raising satellite. This potential entails no angular- 
momentum exchange, and therefore yields no secular effect on the semi-major axis. To prove 
this, let us decompose this potential into a series over the powers of e . This will require of us 
to derive the expression of (a/r)^ . Starting out with the well known development 



a 
r 



= 1 + XI 2-^(^6) cos(pM) 
p=i 

= 1 + e cos(M) + cos(2M) + (^^ cos(3M) - ^ cos(M) ) 

+ (-1 cos(2M) + ^ cos(4M)^ + (^-j^ cos(3M) + ^ cos(M) + ^ cos(5M)^ 



Q 1 "1 "1 \ 

— cos(6M) + — cos(2M) - — cos(4M) j 
one can arrive to the following expansion: 



= l + 6ecos(M)+ + yCos(2M)^e2+ ^^cos(3M) + ^cos(M)^e^ 

/197 101 105\ 4 

+ ( — cos(4M) + — cos(2M) + ^ ) 

/^5157 , 5529 , 2721 , 

+ ( cos(3M) + cos(5M) + cos(M) j 

(483Q 595 739 \ 

^ cos(4M) + 129 cos(2M) + r±l + !£f cos(6M) j + ... , (50) 

whose average over the mean anomaly looks like: 



6 . ^ 15 2 105 4 525 6 , ^, 

: 1 + -e^ + — + — + ... (51) 



Hence, the averaged potential will become: 

In the case of Phobos, the terms of order O(e^) may, in the first approximation, be neglected. 
This means that out of the six Lagrange-type planetary equations the first five will, in the 
first order of e , stay unperturbed, and therefore the elements a , e , a; , ^ , ^2 will, in the first 
order over e , remain unchanged. The Lagrange equation for the longitude will be the only one 
influenced by Uq . That equation will assume the form: 

^ = „ _ A ^ (63) 
at n a da 



which gives 



6 k2 (Mo + m) G m 
L ^ nt + — ^ ° ^ t 54 
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We see that in the first order of e the only secular effect stemming from the potential Uq is a 
linear in time evolution of the longitude. 
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